Section: New Results

Modeling

Local approximate DtN exterior boundary condition

Participants : Hélène Barucq, Rabia Djellouli, Anne-Gaëlle Saint-Guirons.

We investigate analytically the asymptotic behavior of high-order spurious prolate spheroidal modes induced by a second-order local approximate DtN absorbing boundary condition (DtN2) when employed for solving high-frequency acoustic scattering problems. We prove that these reflected modes decay exponentially in the high frequency regime. This theoretical result demonstrates the great potential of the considered absorbing boundary condition for solving efficiently exterior high-frequency Helmholtz problems. In addition, this exponential decay proves the superiority of DtN2 over the widely used Bayliss-Gunsburger-Turkel absorbing boundary condition. This work has been accepted for publication in Progress In Electromagnetics Research B. [19] .

Non-reflecting boundary condition on ellipsoidal boundary

Participants : Hélène Barucq, Anne-Gaëlle Saint-Guirons, Sébastien Tordeux.

The modeling of wave propagation problems using finite element methods usually requires the truncation of the computational domain around the scatterer of interest. Absorbing boundary condition are classically considered in order to avoid spurious reflections. In this paper, we investigate some properties of the Dirichlet to Neumann map posed on a spheroidal boundary in the context of the Helmholtz equation. We focus on the impedance coefficients defining the DtN condition and we aim at establishing suitable properties in order to propose an accurate numerical method for their computation. Then, we state the well-posedness of the corresponding mixed problem and propose a variational formulation adapted to a finite element discretization. This work has been submitted.

A new modified equation approach for solving the wave equation

Participants : Cyril Agut, Hélène Barucq, Julien Diaz, Florent Ventimiglia.

Participants : Roland Martin, Dimitri Komatitsch.

The new method involving $p$-harmonic operator described in section  3.2 has been presented in [13] . We have proved the convergence of the scheme and its stability under a CFL condition. Numerical results in one, two and three-dimensional configurations show that this CFL condition is slightly greater than the CFL condition of the second-order Leap-Frog scheme. We have also studied the penalization parameters involved in the new schemes and their influence of the CFL condition. These results are presented in the PhD. thesis of Cyril Agut [11] . In the framework of the Phd thesis of Florent Ventimiglia, we are now considering the extension of this technique to the first order formulation of the elastodynamic equations.

Stability Analysis of an Interior Penaly Discontinuous Galerkin Method for the Wave equation

Participants : Cyril Agut, Hélène Barucq, Julien Diaz.

The Interior Penalty Discontinuous Galerkin Method  [42] , [38] , [61] we use in the IPDGFEM code requires the introduction of a penalty parameter. Except for regular quadrilateral or cubic meshes, the optimal value of this parameter is not explicitely known. Moreover, the condition number of the resulting stiffness matrix is an increasing function of this parameter, but the precise behaviour has not been explicited neither. We have carried out a theoretical and numerical study of the CFL condition for quadrilateral and cubic meshes, which is presented in the Phd thesis of Cyril Agut [11] . These results were also presented at the peer-reviewed conference Waves 2011 (Vancouver, Canada, July 2011).

Higher Order Absorbing Boundary Conditions for the Wave Equation

Participants : Hélène Barucq, Julien Diaz, Véronique Duprat.

The numerical simulation of wave propagation is generally performed by truncating the propagation medium. Absorbing boundary conditions are then needed. We construct a new family of absorbing boundary conditions from the factorization of the wave equation formulated as a first order system. Using the method of M.E. Taylor, we show that we can generate an infinite number of boundary conditions which can not be obtained via the Niremberg’s factorization method. The conditions can be applied on arbitrarily-shaped surfaces and involve second-order derivatives. We then propose a reduced formulation of the wave equation using an auxiliary unknown which is defined on the regular surface only. The reduced problem allows one to easily include the boundary conditions inside the variational formulation. The corresponding boundary value problem remains well-posed in suitable Hilbert spaces and we give a demonstration in a framework that is suitable to applications. We then study the long-time behavior of the wave field and we show that it tends to 0 as time tends to infinity. This provides a weak stability result that should be completed in the second part of this work. We have then decided to improve the stability result by performing a quantitative study of the energy. We have then shown that the energy is exponentially decaying if the obstacle is star-shaped and the external boundary is convex. This work has been published as INRIA Resarch Reports [34] , [35] and two papers are submitted. We have next addressed the issue of enriching these ABCs by representing evanescent and damping waves. This has given rise to a work for the Hemlholtz equation and we have shown that the enriched ABCs performed better than standard ABCs. The xetension to the acoustic wave has led to new conditions involvin fractional derivatives. To the bast of our knowledge, it is the first time that fractional derivatives have been used for optimizing the performance of ABCs. These new results have been presented in two seminars (University of Bordeaux I and University of Genova) and in two conferences (FACM11 and WAVES 2011). A paper has been published  [17] for the case of evanescent waves and two papers are in preparation. All these results are presented in the PhD thesis of Véronique Duprat [12] .

Numerical methods combining local time stepping and mixed hybrid elements for the terrestrial migration

Participants : Caroline Baldassari, Hélène Barucq, Henri Calandra [Expert Engineer, TOTAL] , Bertrand Denel [Research Engineer, TOTAL] , Julien Diaz, Florent Ventimiglia.

In order to justify the use of our code IPDGFEM for the Reverse Time Migration, we have carried out a performance analysis of the Interior Penalty Discontinuous Galerkin method and of the Spectral Element Method. This analysis, which shows that IPDG performs as well as SEM, has been presented in [14] .

Another aspect of the work concerns the design of local time-stepping algorithms. The local-time stepping strategy proposed in [5] allows for high-order time schemes where the time scheme is adapted to the various space step of the mesh. However, when the mesh contains both low-order and high-order cells, this method not allows for the adaptation of the order of the time-scheme to the order of the cells. We have then presented a new local time-stepping algorithm where both the order of the scheme and the time step vary in the different parts of the mesh. This method has been presented in [14] and at the peer-reviewed conferences Waves 2011 (Vancouver, Canada, July 2011) and DD20 (Domain Decompostion, San Diego, USA, February 2011).

The local-time stepping algorithm is not adapted to handle dissipation terms. A method has been proposed in  [59] , but it is based on an Adams-Bashworth scheme and it requires the storage of additional unknowns. We can not use this scheme for the simulation of seismic waves in very large heterogeneous domains due to memory limitation. We are now working on the design of alternative schemes which would not require the introduction of the auxiliary unknowns. This one of the topics of the PhD. thesis of Florent Ventimiglia.

Perfectly Matched Layers for the Shallow Water equations

Participants : Hélène Barucq, Julien Diaz, Mounir Tlemcani [Assistant Professor, University of Oran, Algeria] .

In  [45] , we have proposed a new Perfectly Matched Layer for Shallow Water equations. This layer required the computation of an auxiliary variable in the whole computational domain. We are now considering a new strategy, which only requires the computation of the auxiliary variable inside the layer. Moeover, the new methodology seems to be well-adapted to the non-linear shallow water equations. We are now performing numerical tests to confirm this point.

Multiperforated plates in linear acoustics

Participants : Abderrahmane Bendali, M'Barek Fares, Sophie Laurens, Estelle Piot, Sébastien Tordeux.

Acoustic engineers use approximate heuristic models to deal with multiperforated plates in liners and in combustion chambers of turbo-engines. These models were suffering from a lack of mathematical justifications and were consequently difficult to improve. Performing an asymptotic analysis (the small parameter is the radius of the perforations), we have justified these models and proposed some improvement. Our theoretical results have been compared to numerical simulations performed at CERFACS (M'Barek Fares) and to acoustical experiments realized at ONERA (Estelle Piot). Two papers are in preparation.

Asymptotic modeling in electromagnetism

Participants : François Buret, Gabriel Caloz, Monique Dauge, Patrick Dular, Marc Duruflé, Erwan Faou, Laurent Krähenbühl, Victor Péron, Ronan Perrussel, Clair Poignard, Damien Voyer.

The following results rely on several problematics developed in section  3.2 , item Asymptotic modeling.

We consider in [21] the equations of electromagnetism set on a domain made of a dielectric and a conductor subdomain in a regime where the conductivity is large. Assuming smoothness for the dielectric–conductor interface, relying on recent works we prove that the solution of the Maxwell equations admits a multiscale asymptotic expansion with profile terms rapidly decaying inside the conductor. This skin effect is measured by introducing a skin depth function that turns out to depend on the mean curvature of the boundary of the conductor. We then confirm these asymptotic results by numerical experiments in various axisymmetric configurations. We also investigate numerically the case of a nonsmooth interface, namely a cylindrical conductor.

We derive new thin layer models in electromagnetism, in [22] . We study the behavior of the electromagnetic field in a biological cell modeled by a medium surrounded by a thin layer and embedded in an ambient medium. We derive approximate transmission conditions in order to replace the membrane by these conditions on the boundary of the interior domain. Our approach is essentially geometric and based on a suitable change of variables in the thin layer. Few notions of differential calculus are given in order to obtain the first-order conditions in a simple way, and numerical simulations validate the theoretical results. Asymptotic transmission conditions at any order are given.

We present a numerical treatment of rounded and sharp corners in the modeling of 2D electrostatic fields in [36] . This work deals with numerical techniques to compute electrostatic fields in devices with rounded corners in 2D situations. The approach leads to the solution of two problems: one on the device where rounded corners are replaced by sharp corners and the other on an unbounded domain representing the shape of the rounded corner after an appropriate rescaling. Both problems are solved using different techniques and numerical results are provided to assess the efficiency and the accuracy of the techniques.

Operator Based Upscaling for Discontinuous Galerkin Methods

Participants : Hélène Barucq, Théophile Chaumont, Julien Diaz, Victor Péron.

Realistic numerical simulations of seimic wave propagation are complicated to handle because they must be performed in strongly heterogeneous media. Two different scales must then be taken into account. Indeed, the medium heterogeneities are very small compared to the characteristic dimensions of the propagation medium. To get acccurate numerical solutions, engineers are then forced to use meshes that match the finest scale representing the heterogeneities. Meshing the whole domain with the fine grid leads then to huge linear systems and the computational cost of the numerical method is then very high. It would be thus very interesting to dispose of a numerical method allowing to represent the heterogenities of the medium accurately while computing on a coarse grid. This is the challenge of multiscale approaches like homoegenization or upscaling. In this work, we use an operator-based upscaling method. Operator-based upscaling methods were first developed for elliptic flow problems (see  [41] ) and then extended to hyperbolic problems (see  [62] , [73] , [72] ). Operator-based upscaling method consists in splitting the solution into a coarse and a fine part. The coarse part is defined on a coarse mesh while the fine part is computed on a fine mesh. In order to speed up calculations, artificial boundary conditions (ABC) are imposed. By enforcing suitable ABCs on the boundary of every cells of the coarse mesh, calculations on the fine grid can be carried out locally. The coarse part is next computed globally on the coarse mesh. Operator-based upscaling methods were so far developed in joint with standard finite element discretisation strategy. In this work, we investigate the idea of combining an operator based upscaling method with discontinous Galerkin finite element methods(DGFEM). To begin with, we have used use the interior penalty method as presented in  [42] for elliptic problems and in  [61] , [60] for the wave equation. This is a quite natural way of addressing this issue because we can use a software package that has been already developed in the team. The first results that we have obtained seem to indicate that an DG operator based uspcaling method could be interesting essentially in case of staionnary problems. Nevertheless, the numerical analysis of the discretized problem must be continued. This work has been initiated during the internship of Theophile Chaumont-Frelet who was a fourth year engineer student at Rouen INSA. A paper dealing with the case of the Laplace operator will be submitted soon.

Discontinuous Galerkin Methods for Seismic Wave Propagation

Participants : Hélène Barucq, Caroline Baldassari, Lionel Boillot, Marie Bonnasse, Julien Diaz, Jérôme Luquel, Vanessa Mattesi, Florent Ventimiglia.

In the framework of our collaboration with Total, we are implementing a Discontinuous Galerkin formulation of the first order elastodynamic wave equations in the plateform Diva which is developed by Total. We consider the formulation proposed in  [54] for isotropic media. During her post-doc, Caroline Baldassi has implemented a three dimensional code with Perfectly Matched Layers for this formulation. Jérôme Luquel has implemented the 2D version of this code during his internship. In the framework of the internship of Marie Bonnasse and the PhD thesis of Lionel Boillot, we have extended the formulation to Vertical Transverse Isotropic and Tilted Transverse Isotropic media in both 2D and 3D. The introduction of Absorbing Boundary Conditions or of PML is still an open problem for these types of media. It is one of the topics of the PhD thesis of Lionel Boillot.

The version of the code that we are using assumed that the properties of the media (density, velocity,...) are constant on each cells of the mesh. Discontinuous Galerkin methods allow for considering more general configurations, where these properties vary as polynomial functions inside each cells. Hence, it is not necessary to define the interfaces between the different media before constructing the mesh. The discontinuities are taken into account directly inside each cells. Moreover, we are able to consider smoothly varying media. In the framework of the internship of Vanessa Mattesi, we have implemented polynomial velocities in a Discontinuous Galerkin formulation. We have compared the results obtained with this method to the one obtained with piecewise constant properties. We have observed that the new formulation was more accurate and that it allowed for a simpler construction of the mesh. However, these gains do not counterbalance the increase of the computational induced by the new method. We have then concluded that considering piecewise constant properties was more appropriate to model seismic wave propagation.

Elastic surface waves in crystals

Participants : José Carcione, Fabio Cavallini, Dimitri Komatitsch, Nathalie Favretto-Cristini.

In [25] , we present a review of wave propagation at the surface of anisotropic media (crystal symmetries). The physics for media of cubic and hexagonal symmetries has been extensively studied based on analytical and semi-analytical methods. However, some controversies regarding surfaces waves and the use of different notations for the same modes require a review of the research done and a clarification of the terminology. In a companion paper we obtain the full-wave solution for the wave propagation at the surface of media with arbitrary symmetry (including cubic and hexagonal symmetries) using two spectral numerical modeling algorithms.

In [27] , we obtain the full-wave solution for the wave propagation at the surface of anisotropic media using two spectral numerical modeling algorithms. The simulations focus on media of cubic and hexagonal symmetries, for which the physics has been reviewed and clarified in a companion paper. Even in the case of homogeneous media, the solution requires the use of numerical methods because the analytical Green's function cannot be obtained in the whole space. The algorithms proposed here allow for a general material variability and the description of arbitrary crystal symmetry at each grid point of the numerical mesh. They are based on high-order spectral approximations of the wave field for computing the spatial derivatives. We test the algorithms by comparison to the analytical solution and obtain the wave field at different faces (stress-free surfaces) of apatite, zinc and copper. Finally, we perform simulations in heterogeneous media, where no analytical solution exists in general, showing that the modeling algorithms can handle large impedance variations at the interface.

Application of an elastoplastic spectral-element method to 3D slope stability analysis

Participants : Hom Nath Gharti, Dimitri Komatitsch, Oye Volker, Roland Martin, Jeroen Tromp.

In [26] , we implement a spectral-element method for 3D time-independent elastoplastic problems in geomechanics. As a first application, we use the method for slope stability analyses ranging from small to large scales. The implementation employs an element-by-element preconditioned conjugate gradient solver for efficient storage. The program accommodates material heterogeneity and complex topography. Either simple or complex water table profiles may be used to assess effects of hydrostatic pressure. Both surface loading and pseudostatic seismic loading are implemented. In order to simulate elastoplastic behavior of slopes, a Mohr-Coulomb yield criterion is employed using an initial strain method (i.e., a viscoplastic algorithm). For large-scale problems, the software is parallelized based on domain decomposition using MPI (Message Passing Interface). Strong-scaling measurements demonstrate that the parallelized software performs efficiently.We validate our spectral-element results against several other methods, and apply the technique to simulate failure of an earthen embankment and a mountain slope.

Indirect Boundary Element Method applied to Fluid-Solid Interfaces

Participants : A. Rodriguez-Castellanos, E. Flores, F. J. Sánchez-Sesma, Carlos Ortiz-Aleman, M. Nava-Flores, Roland Martin.

In [31] , scattering of elastic waves in fluid-solid interfaces is investigated. We use the Indirect Boundary Element Method to study this wave propagation phenomenon in 20 models. Three models are analyzed: a first one with an interface between two half-spaces, one fluid on the top part and the other solid in the bottom; a second model including a fluid half-space above a layered solid; and finally, a third model with a fluid layer bounded by two solid half-spaces. The source, represented by Hankel's function of the second kind, is always applied in the fluid. This indirect formulation can give to the analyst a deep physical insight on the generated diffracted waves because it is closer to the physical reality and can be regarded as a realization of Huygens' principle. In any event, mathematically it is fully equivalent to the classical Somigliana's representation theorem. In order to gauge accuracy we test our method by comparing with an analytical solution known as Discrete Wave Number. A near interface pulse generates scattered waves that can be registered by receivers located in the fluid and it is possible to infer wave velocities of solids. Results are presented in both time and frequency domain, where several aspects related to the different wave types that emerge from this kind of problems are pointed out.

Multiperforated plates in linear acoustics

Participants : Mohamed Amara, Sharang Chaudhry, Julien Diaz, Rabia Djellouli, Magdalena Grigoroscuta-Strugaru.

We have designed a new and efficient solution methodology for solving high-frequency Helmholtz problems. The proposed method is a least-squares based technique that employs variable bases of plane waves at the element level of the domain partition. A local wave tracking strategy is adopted for the selection of the basis at the regional/element level. More specifically, for each element of the mesh partition, a basis of plane waves is chosen so that one of the plane waves in the basis is oriented in the direction of the propagation of the field inside the considered element. The determination of the direction of the field inside the mesh partition is formulated as a minimization problem. Since the problem is nonlinear, we apply Newton's method to determine the minimum. The computation of Jacobians and Hessians that arise in the iterations of the Newtonâs method is based on the exact characterization of the Fréchet derivatives of the field with respect to the propagation directions. Such a characterization is crucial for the stability, fast convergence, and computational efficiency of the Newton algorithm. These results are part of the Master thesis of Sharang Chaudhry (student à CSUN).